module med_map_mod

  use med_kind_mod          , only : CX=>SHR_KIND_CX, CS=>SHR_KIND_CS, CL=>SHR_KIND_CL, R8=>SHR_KIND_R8
  use shr_const_mod         , only : shr_const_pi
  use ESMF                  , only : ESMF_Field, ESMF_FieldGet, ESMF_SUCCESS, ESMF_FAILURE
  use ESMF                  , only : ESMF_LOGMSG_ERROR, ESMF_LOGMSG_INFO, ESMF_LogWrite
  use esmFlds               , only : mapbilnr, mapconsf, mapconsd, mappatch, mapfcopy
  use esmFlds               , only : mapunset, mapnames, nmappers
  use esmFlds               , only : mapnstod, mapnstod_consd, mapnstod_consf
  use esmFlds               , only : ncomps, compatm, compice, compocn, compname
  use esmFlds               , only : mapfcopy, mapconsd, mapconsf, mapnstod
  use esmFlds               , only : mapuv_with_cart3d
  use esmFlds               , only : med_fldList_entry_type
  use esmFlds               , only : fldListFr, fldListTo
  use esmFlds               , only : coupling_mode
  use med_internalstate_mod , only : InternalState
  use med_constants_mod     , only : ispval_mask       => med_constants_ispval_mask
  use med_constants_mod     , only : czero             => med_constants_czero
  use med_constants_mod     , only : dbug_flag         => med_constants_dbug_flag
  use med_utils_mod         , only : chkerr            => med_utils_ChkErr
  use med_utils_mod         , only : memcheck          => med_memcheck
  use med_methods_mod       , only : FB_getFieldN      => med_methods_FB_getFieldN
  use med_methods_mod       , only : FB_init           => med_methods_FB_Init
  use med_methods_mod       , only : FB_reset          => med_methods_FB_Reset
  use med_methods_mod       , only : FB_Clean          => med_methods_FB_Clean
  use med_methods_mod       , only : FB_GetFldPtr      => med_methods_FB_GetFldPtr
  use med_methods_mod       , only : FB_FieldRegrid    => med_methods_FB_FieldRegrid
  use med_methods_mod       , only : FB_Field_diagnose => med_methods_FB_Field_diagnose
  use med_methods_mod       , only : FB_FldChk         => med_methods_FB_FldChk
  use med_methods_mod       , only : FB_GetFieldByName => med_methods_FB_GetFieldByName
  use med_methods_mod       , only : Field_diagnose    => med_methods_Field_diagnose
  use perf_mod              , only : t_startf, t_stopf

  implicit none
  private

  ! public routines
  public :: med_map_RouteHandles_init
  public :: med_map_Fractions_init
  public :: med_map_MapNorm_init
  public :: med_map_FB_Regrid_Norm
  public :: med_map_uv_cart3d

  interface med_map_FB_Regrid_norm
     module procedure med_map_FB_Regrid_Norm_All
     module procedure med_map_FB_Regrid_Norm_Frac
  end interface

  ! private module variables

  character(len=CS)       :: flds_scalar_name
  integer                 :: srcTermProcessing_Value = 0 ! should this be a module variable?
  logical                 :: mastertask
  character(*), parameter :: u_FILE_u  = &
       __FILE__

!================================================================================
contains
!================================================================================

  subroutine med_map_RouteHandles_init(gcomp, llogunit, rc)

    !---------------------------------------------
    ! Initialize route handles in the mediator
    ! Assumptions:
    !   -  Route handles are created per target field bundles NOT
    !      per individual fields in the bundle
    !   -  ALL fields in the bundle are on identical grids
    !   -  MULTIPLE route handles are going to be generated for
    !      given field bundle source and destination grids
    !    - Route handles will ONLY be created if coupling is active
    !      between n1 and n2
    ! Algorithm
    !     n1=source component index
    !     n2=destination component index
    !     nf=field index
    !     fldListFr(n)%flds(nf) is queried to determine the mapindex and mapfile
    !     and the appropriate route handle is created
    !
    ! Regridding is done on a per-field basis AND only for those fields that have a
    ! valid mapping index for the destination component
    !     n = source field index field index
    !     destcomp = destination component index
    !     The fldsSrc input argument is queried for the mapping type of the field
    !     for the desination component
    !        mapindex = fldsSrc(n)%mapindex(destcomp)
    !     If the mapindex is 0 (there is no valid mapping) then NO mapping is done
    !        for the field
    !---------------------------------------------

    use ESMF  , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS, ESMF_LogFlush, ESMF_KIND_I4
    use ESMF  , only : ESMF_GridComp, ESMF_VM, ESMF_Field, ESMF_PoleMethod_Flag, ESMF_POLEMETHOD_ALLAVG
    use ESMF  , only : ESMF_GridCompGet, ESMF_VMGet, ESMF_FieldSMMStore, ESMF_RouteHandleIsCreated
    use ESMF  , only : ESMF_FieldRedistStore, ESMF_FieldRegridStore, ESMF_REGRIDMETHOD_BILINEAR
    use ESMF  , only : ESMF_UNMAPPEDACTION_IGNORE, ESMF_REGRIDMETHOD_CONSERVE, ESMF_NORMTYPE_FRACAREA
    use ESMF  , only : ESMF_REGRIDMETHOD_NEAREST_STOD
    use ESMF  , only : ESMF_NORMTYPE_DSTAREA, ESMF_REGRIDMETHOD_PATCH, ESMF_RouteHandlePrint
    use NUOPC , only : NUOPC_Write

    ! input/output variables
    type(ESMF_GridComp)  :: gcomp
    integer, intent(in)  :: llogunit
    integer, intent(out) :: rc

    ! local variables
    type(InternalState)     :: is_local
    type(ESMF_VM)           :: vm
    type(ESMF_Field)        :: fldsrc
    type(ESMF_Field)        :: flddst
    integer                 :: localPet
    integer                 :: n,n1,n2,m,nf,nflds,ncomp
    integer                 :: SrcMaskValue
    integer                 :: DstMaskValue
    character(len=128)      :: value
    character(len=128)      :: rhname
    character(len=128)      :: rhname_file
    character(len=CS)       :: mapname
    character(len=CX)       :: mapfile
    character(len=CS)       :: string
    integer                 :: mapindex
    logical                 :: rhprint_flag = .false.
    logical                 :: mapexists = .false.
    real(R8)      , pointer :: factorList(:)
    character(CL) , pointer :: fldnames(:)
    !integer(ESMF_KIND_I4), pointer :: unmappedDstList(:)
    character(len=128)      :: logMsg
    type(ESMF_PoleMethod_Flag), parameter :: polemethod=ESMF_POLEMETHOD_ALLAVG
    character(len=*), parameter :: subname=' (module_med_map: RouteHandles_init) '
    !-----------------------------------------------------------
    call t_startf('MED:'//subname)

    if (dbug_flag > 1) then
       call ESMF_LogWrite("Starting to initialize RHs", ESMF_LOGMSG_INFO)
       call ESMF_LogFlush()
    endif

    rc = ESMF_SUCCESS

    ! Determine mastertask
    call ESMF_GridCompGet(gcomp, vm=vm, rc=rc)
    call ESMF_VMGet(vm, localPet=localPet, rc=rc)
    mastertask = .false.
    if (localPet == 0) mastertask=.true.
    ! Get the internal state from Component.
    nullify(is_local%wrap)
    call ESMF_GridCompGetInternalState(gcomp, is_local, rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! Create the necessary route handles
    if (mastertask) write(llogunit,*) ' '
    do n1 = 1, ncomps
       do n2 = 1, ncomps

          if (trim(coupling_mode) == 'cesm') then
             dstMaskValue = ispval_mask
             srcMaskValue = ispval_mask
             if (n1 == compocn .or. n1 == compice) srcMaskValue = 0
             if (n2 == compocn .or. n2 == compice) dstMaskValue = 0
          else if (coupling_mode(1:5) == 'nems_') then
             if (n1 == compatm .and. (n2 == compocn .or. n2 == compice)) then
                srcMaskValue = 1
                dstMaskValue = 0
             else if (n2 == compatm .and. (n1 == compocn .or. n1 == compice)) then
                srcMaskValue = 0
                dstMaskValue = 1
             else if ((n1 == compocn .and. n2 == compice) .or. (n1 == compice .and. n2 == compocn)) then
                srcMaskValue = 0
                dstMaskValue = 0
             else
                ! TODO: what should the condition be here?
                dstMaskValue = ispval_mask
                srcMaskValue = ispval_mask
             end if
          end if

          !--- get single fields from bundles
          !--- 1) ASSUMES all fields in the bundle are on identical grids
          !--- 2) MULTIPLE route handles are going to be generated for
          !---    given field bundle source and destination grids

          if (n1 /= n2) then

             ! Determine route handle names
             rhname = trim(compname(n1))//"2"//trim(compname(n2))

             if (is_local%wrap%med_coupling_active(n1,n2)) then ! If coupling is active between n1 and n2

                call FB_GetFieldN(is_local%wrap%FBImp(n1,n1), 1, fldsrc, rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return

                call FB_GetFieldN(is_local%wrap%FBImp(n1,n2), 1, flddst, rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return

                ! Loop over fields
                do nf = 1,size(fldListFr(n1)%flds)

                   ! Determine the mapping type for mapping field nf from n1 to n2
                   mapindex = fldListFr(n1)%flds(nf)%mapindex(n2)

                   ! separate check first since Fortran does not have short-circuit evaluation
                   if (mapindex == mapunset) cycle

                   ! Create route handle for target mapindex if route handle is required
                   ! (i.e. mapindex /= mapunset) and route handle has not already been created
                   mapexists = .false.
                   if (mapindex == mapnstod_consd .and. &
                        ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapnstod), rc=rc) .and. &
                        ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapconsd), rc=rc)) then
                      mapexists = .true.
                   else if (mapindex == mapnstod_consf .and. &
                        ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapnstod), rc=rc) .and. &
                        ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapconsf), rc=rc)) then
                      mapexists = .true.
                   else if (ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapindex), rc=rc)) then
                      mapexists = .true.
                   end if
                   if (chkerr(rc,__LINE__,u_FILE_u)) return

                   if (.not. mapexists) then

                      mapname  = trim(mapnames(mapindex))
                      mapfile  = trim(fldListFr(n1)%flds(nf)%mapfile(n2))
                      string   = trim(rhname)//'_weights'

                      if (mapindex == mapfcopy) then
                         if (mastertask) then
                            write(llogunit,'(3A)') subname,trim(string),' RH redist '
                         end if
                         call ESMF_LogWrite(trim(subname) // trim(string) // ' RH redist ', ESMF_LOGMSG_INFO)
                         call ESMF_FieldRedistStore(fldsrc, flddst, &
                              routehandle=is_local%wrap%RH(n1,n2,mapindex), &
                              ignoreUnmatchedIndices = .true., rc=rc)
                         if (chkerr(rc,__LINE__,u_FILE_u)) return
                      else if (mapfile /= 'unset') then
                         if (mastertask) then
                            write(llogunit,'(4A)') subname,trim(string),' RH '//trim(mapname)//' via input file ',&
                                 trim(mapfile)
                         end if
                         call ESMF_LogWrite(subname // trim(string) //&
                              ' RH '//trim(mapname)//' via input file '//trim(mapfile), ESMF_LOGMSG_INFO)
                         call ESMF_FieldSMMStore(fldsrc, flddst, mapfile, &
                              routehandle=is_local%wrap%RH(n1,n2,mapindex), &
                              ignoreUnmatchedIndices=.true., &
                              srcTermProcessing=srcTermProcessing_Value, rc=rc)
                         if (chkerr(rc,__LINE__,u_FILE_u)) return
                      else
                         ! Create route handle on the fly
                         if (mastertask) write(llogunit,'(3A)') subname,trim(string),&
                              ' RH regrid for '//trim(mapname)//' computed on the fly'
                         call ESMF_LogWrite(subname // trim(string) //&
                              ' RH regrid for '//trim(mapname)//' computed on the fly', ESMF_LOGMSG_INFO)
                         if (mapindex == mapbilnr) then
                            srcTermProcessing_Value = 0
                            call ESMF_FieldRegridStore(fldsrc, flddst, &
                                 routehandle=is_local%wrap%RH(n1,n2,mapindex), &
                                 srcMaskValues=(/srcMaskValue/), &
                                 dstMaskValues=(/dstMaskValue/), &
                                 regridmethod=ESMF_REGRIDMETHOD_BILINEAR, &
                                 polemethod=polemethod, &
                                 srcTermProcessing=srcTermProcessing_Value, &
                                 factorList=factorList, &
                                 ignoreDegenerate=.true., &
                                 unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
                         else if ((mapindex == mapconsf .or. mapindex == mapnstod_consf) .and. &
                              .not. ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapconsf))) then
                            call ESMF_FieldRegridStore(fldsrc, flddst, &
                                 routehandle=is_local%wrap%RH(n1,n2,mapconsf), &
                                 srcMaskValues=(/srcMaskValue/), &
                                 dstMaskValues=(/dstMaskValue/), &
                                 regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
                                 normType=ESMF_NORMTYPE_FRACAREA, &
                                 srcTermProcessing=srcTermProcessing_Value, &
                                 factorList=factorList, &
                                 ignoreDegenerate=.true., &
                                 unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
                                 !unmappedDstList=unmappedDstList, &
                                 rc=rc)
                         else if ((mapindex == mapconsd .or. mapindex == mapnstod_consd) .and. &
                            .not. ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapconsd))) then
                            call ESMF_FieldRegridStore(fldsrc, flddst, &
                                 routehandle=is_local%wrap%RH(n1,n2,mapconsd), &
                                 srcMaskValues=(/srcMaskValue/), &
                                 dstMaskValues=(/dstMaskValue/), &
                                 regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
                                 normType=ESMF_NORMTYPE_DSTAREA, &
                                 srcTermProcessing=srcTermProcessing_Value, &
                                 factorList=factorList, &
                                 ignoreDegenerate=.true., &
                                 unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
                                 !unmappedDstList=unmappedDstList, &
                                 rc=rc)
                         else if (mapindex == mappatch) then
                            call ESMF_FieldRegridStore(fldsrc, flddst, &
                                 routehandle=is_local%wrap%RH(n1,n2,mapindex), &
                                 srcMaskValues=(/srcMaskValue/), &
                                 dstMaskValues=(/dstMaskValue/), &
                                 regridmethod=ESMF_REGRIDMETHOD_PATCH, &
                                 polemethod=polemethod, &
                                 srcTermProcessing=srcTermProcessing_Value, &
                                 factorList=factorList, &
                                 ignoreDegenerate=.true., &
                                 unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
                         end if
                         ! consd_nstod method requires a second routehandle
                         if ((mapindex == mapnstod .or. mapindex == mapnstod_consd .or. mapindex == mapnstod_consf) .and. &
                              .not. ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapnstod),rc=rc)) then
                            call ESMF_FieldRegridStore(fldsrc, flddst, &
                                 routehandle=is_local%wrap%RH(n1,n2,mapnstod), &
                                 srcMaskValues=(/srcMaskValue/), &
                                 dstMaskValues=(/dstMaskValue/), &
                                 regridmethod=ESMF_REGRIDMETHOD_NEAREST_STOD, &
                                 srcTermProcessing=srcTermProcessing_Value, &
                                 factorList=factorList, &
                                 ignoreDegenerate=.true., &
                                 unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, &
                                 rc=rc)
                         end if
                         if (chkerr(rc,__LINE__,u_FILE_u)) return
                         if (rhprint_flag .and. mapindex /= mapnstod_consd .and. mapindex /= mapnstod_consf) then
                            call NUOPC_Write(factorList, "array_med_"//trim(string)//"_consf.nc", rc)
                            if (chkerr(rc,__LINE__,u_FILE_u)) return
                         end if
                         !if (associated(unmappedDstList)) then
                         !   write(logMsg,*) trim(subname),trim(string),'     number of unmapped dest points = ', size(unmappedDstList)
                         !   call ESMF_LogWrite(trim(logMsg), ESMF_LOGMSG_INFO)
                         !end if
                      end if
                      if (rhprint_flag .and. mapindex /= mapnstod_consd .and. mapindex /= mapnstod_consf) then
                         call ESMF_LogWrite(trim(subname)//trim(string)//": printing  RH for "//trim(mapname), &
                              ESMF_LOGMSG_INFO)
                         call ESMF_RouteHandlePrint(is_local%wrap%RH(n1,n2,mapindex), rc=rc)
                         if (chkerr(rc,__LINE__,u_FILE_u)) return
                      endif
                      if (chkerr(rc,__LINE__,u_FILE_u)) return
                      ! Check that a valid route handle has been created
                      if (  mapindex /= mapnstod_consd .and. mapindex /= mapnstod_consf .and. &
                           .not. ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,mapindex), rc=rc)) then
                         call ESMF_LogWrite(trim(subname)//trim(string)//": failed   RH "//trim(mapname), &
                              ESMF_LOGMSG_INFO)
                      endif
                   end if
                end do ! loop over fields
             end if ! if coupling is active between n1 and n2
          end if ! if n1 not equal to n2
       end do ! loop over n2
    end do ! loop over n1

    if (dbug_flag > 1) then
      call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO)
    endif
    call t_stopf('MED:'//subname)

  end subroutine med_map_RouteHandles_init

  !================================================================================

  subroutine med_map_Fractions_init(gcomp, n1, n2, FBSrc, FBDst, RouteHandle, rc)

    !---------------------------------------------
    ! Initialize initialize additional route handles
    ! for mapping fractions
    !---------------------------------------------

    use ESMF  , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS, ESMF_LogFlush
    use ESMF  , only : ESMF_GridComp, ESMF_FieldBundle, ESMF_RouteHandle, ESMF_Field
    use ESMF  , only : ESMF_FieldRedistStore, ESMF_FieldSMMStore, ESMF_FieldRegridStore
    use ESMF  , only : ESMF_UNMAPPEDACTION_IGNORE, ESMF_REGRIDMETHOD_CONSERVE, ESMF_NORMTYPE_FRACAREA
    use NUOPC , only : NUOPC_CompAttributeGet

    type(ESMF_GridComp)                    :: gcomp
    integer                , intent(in)    :: n1
    integer                , intent(in)    :: n2
    type(ESMF_FieldBundle) , intent(in)    :: FBSrc
    type(ESMF_FieldBundle) , intent(in)    :: FBDst
    type(ESMF_RouteHandle) , intent(inout) :: RouteHandle
    integer                , intent(out)   :: rc

    ! local variables
    type(ESMF_Field)   :: fldsrc
    type(ESMF_Field)   :: flddst
    character(len=128) :: rhname
    character(len=CS)  :: mapname
    character(len=CX)  :: mapfile
    character(len=CS)  :: string
    integer            :: SrcMaskValue
    integer            :: DstMaskValue
    real(R8), pointer  :: factorList(:)
    character(len=*), parameter :: subname=' (med_map_fractions_init: ) '
    !---------------------------------------------
    call t_startf('MED:'//subname)

    if (dbug_flag > 1) then
       call ESMF_LogWrite("Initializing RHs not yet created and needed for mapping fractions", &
            ESMF_LOGMSG_INFO)
       call ESMF_LogFlush()
    endif

    rc = ESMF_SUCCESS

    call FB_getFieldN(FBsrc, 1, fldsrc, rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    call FB_getFieldN(FBDst, 1, flddst, rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    dstMaskValue = ispval_mask
    srcMaskValue = ispval_mask
    if (n1 == compocn .or. n1 == compice) srcMaskValue = 0
    if (n2 == compocn .or. n2 == compice) dstMaskValue = 0

    rhname = trim(compname(n1))//"2"//trim(compname(n2))
    string   = trim(rhname)//'_weights'
    if ( (n1 == compocn .and. n2 == compice) .or. (n1 == compice .and. n2 == compocn)) then
       mapfile = 'idmap'
    else
       call ESMF_LogWrite("Querying for attribute "//trim(rhname)//"_fmapname = ", ESMF_LOGMSG_INFO)
       call NUOPC_CompAttributeGet(gcomp, name=trim(rhname)//"_fmapname", value=mapfile, rc=rc)
       mapname = trim(mapnames(mapconsf))
    end if

    if (mapfile == 'idmap') then
       call ESMF_LogWrite(trim(subname) // trim(string) //&
            ' RH '//trim(mapname)// ' is redist', ESMF_LOGMSG_INFO)
       call ESMF_FieldRedistStore(fldsrc, flddst, &
            routehandle=RouteHandle, &
            ignoreUnmatchedIndices = .true., rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    else if (mapfile /= 'unset') then
       call ESMF_LogWrite(subname // trim(string) //&
            ' RH '//trim(mapname)//' via input file '//trim(mapfile), ESMF_LOGMSG_INFO)
       call ESMF_FieldSMMStore(fldsrc, flddst, mapfile, &
            routehandle=RouteHandle, &
            ignoreUnmatchedIndices=.true., &
            srcTermProcessing=srcTermProcessing_Value, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    else
       call ESMF_LogWrite(subname // trim(string) //&
            ' RH '//trim(mapname)//' computed on the fly '//trim(mapfile), ESMF_LOGMSG_INFO)
       call ESMF_FieldRegridStore(fldsrc, flddst, &
            routehandle=RouteHandle, &
            srcMaskValues=(/srcMaskValue/), &
            dstMaskValues=(/dstMaskValue/), &
            regridmethod=ESMF_REGRIDMETHOD_CONSERVE, &
            normType=ESMF_NORMTYPE_FRACAREA, &
            srcTermProcessing=srcTermProcessing_Value, &
            factorList=factorList, &
            ignoreDegenerate=.true., &
            unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)
    end if

    if (dbug_flag > 1) then
      call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO)
    endif
    call t_stopf('MED:'//subname)

  end subroutine med_map_Fractions_init

!================================================================================

  subroutine med_map_MapNorm_init(gcomp, llogunit, rc)

    !---------------------------------------
    ! Initialize unity normalization field bundle
    ! and do the mapping for unity normalization up front
    !---------------------------------------

    use ESMF , only: ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS, ESMF_LogFlush
    use ESMF , only: ESMF_GridComp, ESMF_FieldBundle, ESMF_RouteHandleIsCreated

    ! input/output variables
    type(ESMF_GridComp)  :: gcomp
    integer, intent(in)  :: llogunit
    integer, intent(out) :: rc

    ! local variables
    type(InternalState)        :: is_local
    type(ESMF_FieldBundle)     :: FBTmp
    integer                    :: n1, n2, m
    character(len=CS)          :: normname
    character(len=1)           :: cn1,cn2,cm
    real(R8), pointer          :: dataptr(:)
    character(len=*),parameter :: subname='(module_MED_MAP:MapNorm_init)'
    !-----------------------------------------------------------
    call t_startf('MED:'//subname)

    if (dbug_flag > 1) then
       call ESMF_LogWrite("Starting to initialize unity map normalizations", ESMF_LOGMSG_INFO)
       call ESMF_LogFlush()
    endif

    rc = ESMF_SUCCESS

    ! Get the internal state from Component.
    nullify(is_local%wrap)
    call ESMF_GridCompGetInternalState(gcomp, is_local, rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! Initialize module variables
    flds_scalar_name = is_local%wrap%flds_scalar_name

    ! Create the normalization field bundles
    normname = 'one'
    do n1 = 1,ncomps
       do n2 = 1,ncomps
          if (n1 /= n2) then
             if (is_local%wrap%med_coupling_active(n1,n2)) then
                do m = 1,nmappers
                   if (ESMF_RouteHandleIsCreated(is_local%wrap%RH(n1,n2,m), rc=rc)) then
                      if (dbug_flag > 1) then
                         write(cn1,'(i1)') n1; write(cn2,'(i1)') n2; write(cm ,'(i1)') m
                         call ESMF_LogWrite(trim(subname)//":"//'creating FBMapNormOne for '&
                              //compname(n1)//'->'//compname(n2)//'with mapping '//mapnames(m), &
                              ESMF_LOGMSG_INFO)
                      endif
                      call FB_init(FBout=is_local%wrap%FBNormOne(n1,n2,m), &
                           flds_scalar_name=flds_scalar_name, &
                           FBgeom=is_local%wrap%FBImp(n1,n2), &
                           fieldNameList=(/trim(normname)/), name='FBNormOne', rc=rc)
                      if (chkerr(rc,__LINE__,u_FILE_u)) return

                      call FB_reset(is_local%wrap%FBNormOne(n1,n2,m), value=czero, rc=rc)
                      if (chkerr(rc,__LINE__,u_FILE_u)) return

                      call FB_init(FBout=FBTmp, &
                           flds_scalar_name=flds_scalar_name, &
                           STgeom=is_local%wrap%NStateImp(n1), &
                           fieldNameList=(/trim(normname)/), name='FBTmp', rc=rc)
                      if (chkerr(rc,__LINE__,u_FILE_u)) return

                      call FB_GetFldPtr(FBTmp, trim(normname), fldptr1=dataPtr, rc=rc)
                      if (chkerr(rc,__LINE__,u_FILE_u)) return
                      dataptr(:) = 1.0_R8

                      call FB_FieldRegrid(&
                           FBTmp                           , normname, &
                           is_local%wrap%FBNormOne(n1,n2,m), normname, &
                           is_local%wrap%RH(n1,n2,m), rc)
                      if (chkerr(rc,__LINE__,u_FILE_u)) return

                      call FB_clean(FBTmp, rc=rc)
                      if (chkerr(rc,__LINE__,u_FILE_u)) return
                   end if
                end do
             end if
          end if
       end do
    end do

    if (dbug_flag > 1) then
      call ESMF_LogWrite(trim(subname)//": done", ESMF_LOGMSG_INFO)
    endif
    call t_stopf('MED:'//subname)

  end subroutine med_map_MapNorm_init

  !================================================================================

  subroutine med_map_FB_Regrid_Norm_All(fldsSrc, srccomp, destcomp, &
       FBSrc, FBDst, FBFracSrc, FBNormOne, RouteHandles, string, rc)

    ! ----------------------------------------------
    ! Map field bundles with appropriate fraction weighting
    ! ----------------------------------------------

    use NUOPC    , only: NUOPC_IsConnected
    use ESMF     , only: ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS
    use ESMF     , only: ESMF_LOGMSG_ERROR, ESMF_FAILURE, ESMF_MAXSTR
    use ESMF     , only: ESMF_Mesh, ESMF_MeshGet, ESMF_MESHLOC_ELEMENT, ESMF_TYPEKIND_R8
    use ESMF     , only: ESMF_FieldBundle, ESMF_FieldBundleIsCreated, ESMF_FieldBundleGet
    use ESMF     , only: ESMF_RouteHandle, ESMF_RouteHandleIsCreated
    use ESMF     , only: ESMF_REGION_SELECT, ESMF_REGION_TOTAL
    use ESMF     , only: ESMF_Field, ESMF_FieldGet, ESMF_FieldIsCreated
    use ESMF     , only: ESMF_FieldDestroy, ESMF_FieldCreate
    use ESMF     , only: ESMF_TERMORDER_SRCSEQ, ESMF_Region_Flag, ESMF_REGION_TOTAL
    use ESMF     , only: ESMF_REGION_SELECT, ESMF_FieldRegrid

    ! input/output variables
    type(med_fldList_entry_type) , pointer       :: fldsSrc(:)
    integer                            , intent(in)    :: srccomp
    integer                            , intent(in)    :: destcomp
    type(ESMF_FieldBundle)             , intent(inout) :: FBSrc
    type(ESMF_FieldBundle)             , intent(inout) :: FBDst
    type(ESMF_FieldBundle)             , intent(in)    :: FBFracSrc
    type(ESMF_FieldBundle)             , intent(in)    :: FBNormOne(:)
    type(ESMF_RouteHandle)             , intent(inout) :: RouteHandles(:)
    character(len=*), optional         , intent(in)    :: string
    integer                            , intent(out)   :: rc

    ! local variables
    integer               :: i, n, k
    integer               :: lrank
    character(len=CS)     :: lstring
    integer               :: mapindex
    character(len=CS)     :: mapnorm
    character(len=CS)     :: fldname
    type(ESMF_Mesh)       :: lmesh
    type(ESMF_Field)      :: srcField
    type(ESMF_Field)      :: dstField
    type(ESMF_Field)      :: lfield
    type(ESMF_Field)      :: frac_field_src
    type(ESMF_Field)      :: frac_field_dst
    real(R8), allocatable :: data_srctmp(:)
    real(R8), allocatable :: data_srctmp_1d(:)
    real(R8), allocatable :: data_srctmp_2d(:,:)  
    real(R8), pointer     :: data_src_1d(:)
    real(R8), pointer     :: data_src_2d(:,:)
    real(R8), pointer     :: data_frac(:)
    real(R8), pointer     :: data_norm(:)
    logical               :: used_cart3d_for_uvmapping
    logical               :: frac_field_created
    type(ESMF_Field)      :: usrc,vsrc
    type(ESMF_Field)      :: udst,vdst
    integer               :: ungriddedUBound(1)     ! currently the size must equal 1 for rank 2 fields
    integer               :: gridToFieldMap(1)      ! currently the size must equal 1 for rank 2 fields
    logical               :: checkflag = .false.
    character(len=*), parameter :: subname='(module_MED_Map:med_map_Regrid_Norm)'
    !-------------------------------------------------------------------------------

    call t_startf('MED:'//subname)
    call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO)
    call memcheck(subname, 1, mastertask)

#ifdef DEBUG
    checkflag = .true.
#endif

    !---------------------------------------

    if (present(string)) then
      lstring = trim(string)
    else
      lstring = " "
    endif

    rc = ESMF_SUCCESS

    !---------------------------------------
    ! First - reset the field bundle on the destination grid to zero
    !---------------------------------------

    call FB_reset(FBDst, value=czero, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    !---------------------------------------
    ! Loop over all fields in the source field bundle and map them to
    ! the destination field bundle accordingly
    !---------------------------------------

    call ESMF_LogWrite(trim(subname)//" *** mapping from "//trim(compname(srccomp))//" to "//&
         trim(compname(destcomp))//" ***", ESMF_LOGMSG_INFO)

    frac_field_created = .false.
    used_cart3d_for_uvmapping = .false.
    do n = 1,size(fldsSrc)
       ! Determine if field is a scalar - and if so go to next iternation
       fldname  = fldsSrc(n)%shortname
       if (fldname == flds_scalar_name) CYCLE

       ! Determine if there is a map index and if its zero go to next iteration
       mapindex = fldsSrc(n)%mapindex(destcomp)
       if (mapindex == 0) CYCLE
       mapnorm  = fldsSrc(n)%mapnorm(destcomp)

       ! Determine if field is FBSrc or FBDst or connected - and if not go to next iteration
       if (.not. FB_FldChk(FBSrc, trim(fldname), rc=rc)) then
          if (dbug_flag > 5) then
             call ESMF_LogWrite(trim(subname)//" field not found in FBSrc: "//trim(fldname), ESMF_LOGMSG_INFO)
          end if
          CYCLE
       else if (.not. FB_FldChk(FBDst, trim(fldname), rc=rc)) then
          if (dbug_flag > 5) then
             call ESMF_LogWrite(trim(subname)//" field not found in FBDst: "//trim(fldname), ESMF_LOGMSG_INFO)
          end if
          CYCLE
       end if

       ! -------------------
       ! Error checks
       ! -------------------

       if (.not. FB_FldChk(FBSrc, fldname, rc=rc)) then
          call ESMF_LogWrite(trim(subname)//" field not found in FBSrc: "//trim(fldname), ESMF_LOGMSG_INFO)
       else if (.not. FB_FldChk(FBDst, fldname, rc=rc)) then
          call ESMF_LogWrite(trim(subname)//" field not found in FBDst: "//trim(fldname), ESMF_LOGMSG_INFO)
       else if (mapindex == mapnstod_consd) then
          if (.not. ESMF_RouteHandleIsCreated(RouteHandles(mapconsd), rc=rc) .or. &
               .not. ESMF_RouteHandleIsCreated(RouteHandles(mapnstod), rc=rc)) then
             call ESMF_LogWrite(trim(subname)//trim(lstring)//&
                  ": ERROR RH not available for "//mapnames(mapindex)//": fld="//trim(fldname), &
                  ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
             rc = ESMF_FAILURE
             return
          end if
       else if (mapindex == mapnstod_consf) then
          if (.not. ESMF_RouteHandleIsCreated(RouteHandles(mapconsf), rc=rc) .or. &
              .not. ESMF_RouteHandleIsCreated(RouteHandles(mapnstod), rc=rc)) then
             call ESMF_LogWrite(trim(subname)//trim(lstring)//&
                  ": ERROR RH not available for "//mapnames(mapindex)//": fld="//trim(fldname), &
                  ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
             rc = ESMF_FAILURE
             return
          end if
       else if (.not. ESMF_RouteHandleIsCreated(RouteHandles(mapindex), rc=rc)) then
          call ESMF_LogWrite(trim(subname)//trim(lstring)//&
               ": ERROR RH not available for "//mapnames(mapindex)//": fld="//trim(fldname), &
               ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
          rc = ESMF_FAILURE
          return
       end if

       ! -------------------
       ! Do cart3d mapping for u and v fields from atm if appropriate
       ! -------------------

       if (mapuv_with_cart3d) then
          if ((trim(fldname) == 'Sa_u' .or. trim(fldname) == 'Sa_v')) then
             if  (.not. used_cart3d_for_uvmapping) then
                mapindex = fldsSrc(n)%mapindex(destcomp)
                mapnorm  = fldsSrc(n)%mapnorm(destcomp)
                call ESMF_FieldBundleGet(FBSrc, fieldName='Sa_u', field=usrc, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
                call ESMF_FieldBundleGet(FBSrc, fieldName='Sa_v', field=vsrc, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
                call ESMF_FieldBundleGet(FBDst, fieldName='Sa_u', field=udst, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
                call ESMF_FieldBundleGet(FBDst, fieldName='Sa_v', field=vdst, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return

                call ESMF_LogWrite(trim(subname)//" --> remapping "//trim(fldname)//" with "//trim(mapnames(mapindex)), &
                     ESMF_LOGMSG_INFO)

                call med_map_uv_cart3d(usrc, vsrc, udst, vdst, RouteHandles(mapindex), rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return

                used_cart3d_for_uvmapping = .true.
             end if
             CYCLE
          end if
       end if

       ! -------------------
       ! Get the source and destination fields
       ! -------------------

       call ESMF_LogWrite(trim(subname)//" --> remapping "//trim(fldname)//" with "//trim(mapnames(mapindex)), &
            ESMF_LOGMSG_INFO)

       call ESMF_FieldBundleGet(FBSrc, fieldName=trim(fldname), field=srcfield, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       call ESMF_FieldBundleGet(FBDst, fieldName=trim(fldname), field=dstfield, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       ! -------------------
       ! Do the mapping
       ! -------------------

       if (mapindex == mapfcopy) then
          call FB_FieldRegrid(FBSrc, fldname, FBDst, fldname, RouteHandles(mapindex), rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return

       else
          ! Determine the normalization for the map
          mapnorm  = fldsSrc(n)%mapnorm(destcomp)

          if ( trim(mapnorm) /= 'unset' .and. trim(mapnorm) /= 'one' .and. trim(mapnorm) /= 'none') then

             call FB_getFieldByName(FBSrc, fldname, lfield, rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return
             call ESMF_FieldGet(lfield, rank=lrank, rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return

             ! get a pointer to source field data in FBSrc
             if (lrank == 1) then
                call ESMF_FieldGet(srcfield, farrayPtr=data_src_1d, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
             else if (lrank == 2) then
                call ESMF_FieldGet(srcfield, ungriddedUBound=ungriddedUBound, gridToFieldMap=gridToFieldMap, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
                if (gridToFieldMap(1) /= 2) then
                   call ESMF_LogWrite(trim(subname)//" fldname= "//trim(fldname)//&
                        "has gridTofieldMap not equal to 2",  ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
                   rc = ESMF_FAILURE
                   return
                end if
                call ESMF_FieldGet(srcfield, farrayPtr=data_src_2d, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
             end if

             ! allocate memory for a save array if not already allocated
             if (lrank == 1) then
                if (.not. allocated(data_srctmp_1d) .or. size(data_srctmp_1d) /= size(data_src_1d)) then
                   if (allocated(data_srctmp_1d)) then
                      deallocate(data_srctmp_1d)
                   endif
                   allocate(data_srctmp_1d(size(data_src_1d)))
                endif
             elseif (lrank == 2) then
                if (.not. allocated(data_srctmp_2d) .or. size(data_srctmp_2d) /= size(data_src_2d)) then
                   if (allocated(data_srctmp_2d)) then
                      deallocate(data_srctmp_2d)
                   endif
                   allocate(data_srctmp_2d(size(data_src_2d,dim=1), size(data_src_2d,dim=2)))
                endif
             end if

             ! get a pointer to the array of the normalization on the source grid - this must
             ! be the same size is as fraction on the source grid
             call FB_GetFldPtr(FBFracSrc, trim(mapnorm), data_frac, rc=rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return

             ! copy data_src to data_srctmp then multiply data_src by fraction
             if (lrank == 1) then
                data_srctmp_1d(:) = data_src_1d(:)
                data_src_1d(:) = data_src_1d(:) * data_frac(:)
             elseif (lrank == 2) then
                if (size(data_frac) /= size(data_src_2d,dim=2)) then
                   write(6,*)'ERROR: size(frac)   = ',size(data_frac)
                   write(6,*)'ERROR: size(data_src_2d,dim=1),size(data_src_2d,dim=2) = ',&
                        size(data_src_2d,dim=1),size(data_src_2d,dim=2)
                   call ESMF_LogWrite(trim(subname)//" fldname= "//trim(fldname)//&
                        "size of frac not equal to size of distributed data", &
                        ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
                   rc = ESMF_FAILURE
                   return
                end if
                data_srctmp_2d(:,:) = data_src_2d(:,:)
                do i = 1,size(data_frac)
                   data_src_2d(:,i) =  data_src_2d(:,i) * data_frac(i)
                end do
             end if

             ! regrid field with name fldname from FBsrc to FBDst
             call map_field_src2dst (fldname, srcfield, dstfield, RouteHandles, mapindex, rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return

             ! restore original value
             if (lrank == 1) then
                data_src_1d(:) = data_srctmp_1d(:)
             elseif (lrank == 2) then
                data_src_2d(:,:) = data_srctmp_2d(:,:) 
             end if

             ! regrid fraction from source to dest
             if (.not. ESMF_FieldIsCreated(frac_field_dst)) then
                ! get fraction field on source mesh
                call ESMF_FieldBundleGet(FBFracSrc, mapnorm, field=frac_field_src, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return

                ! create fraction field on destination mesh
                call ESMF_FieldBundleGet(FBDst, fldname, field=lfield, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
                call ESMF_FieldGet(lfield, mesh=lmesh, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
                frac_field_dst = ESMF_FieldCreate(lmesh, ESMF_TYPEKIND_R8, name=mapnorm, meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return

                ! regrid fraction field from source to destination
                if (mapindex == mapnstod_consd) then
                   call ESMF_FieldRegrid(frac_field_src, frac_field_dst, RouteHandles(mapnstod), &
                        termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
                   if (chkerr(rc,__LINE__,u_FILE_u)) return
                   call ESMF_FieldRegrid(frac_field_src, frac_field_dst, RouteHandles(mapconsd), &
                        termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_SELECT, rc=rc)
                   if (chkerr(rc,__LINE__,u_FILE_u)) return
                else if (mapindex == mapnstod_consf) then
                   call ESMF_FieldRegrid(frac_field_src, frac_field_dst, RouteHandles(mapnstod), &
                        termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
                   if (chkerr(rc,__LINE__,u_FILE_u)) return
                   call ESMF_FieldRegrid(frac_field_src, frac_field_dst, RouteHandles(mapconsf), &
                        termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_SELECT, rc=rc)
                   if (chkerr(rc,__LINE__,u_FILE_u)) return
                else
                   call ESMF_FieldRegrid(frac_field_src, frac_field_dst, routehandle=RouteHandles(mapindex), &
                        termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
                   if (chkerr(rc,__LINE__,u_FILE_u)) return
                end if

                ! get pointer to mapped fraction
                call ESMF_FieldGet(frac_field_dst, farrayPtr=data_norm, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
             end if

             ! normalize destination mapped values by the reciprocal of the mapped fraction
             call norm_field_dest(trim(fldname), dstfield, data_norm, rc)

             if (dbug_flag > 1) then
                call FB_Field_diagnose(FBDst, fldname, " --> after frac: ", rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
             end if

          else if (trim(mapnorm) == 'one' .or. trim(mapnorm) == 'none') then

             !-------------------------------------------------
             ! unity or no normalization
             !-------------------------------------------------

             ! map source field to destination grid
             call map_field_src2dst (trim(fldname), srcfield, dstfield, RouteHandles, mapindex, rc)
             if (chkerr(rc,__LINE__,u_FILE_u)) return

             ! obtain unity normalization factor and multiply interpolated field by reciprocal of normalization factor
             if (trim(mapnorm) == 'one') then
                call ESMF_FieldBundleGet(FBNormOne(mapindex), fieldName='one', field=lfield, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return
                call ESMF_FieldGet(lfield, farrayPtr=data_norm, rc=rc)
                if (chkerr(rc,__LINE__,u_FILE_u)) return

                call norm_field_dest(trim(fldname), dstfield, data_norm, rc)
             end if ! mapnorm is 'one'

          end if ! mapnorm is 'one' or 'nne'
       end if ! mapindex is not mapfcopy and field exists

       if (dbug_flag > 1) then
          call FB_Field_diagnose(FBDst, fldname, &
               string=trim(subname) //' FBImp('//trim(compname(srccomp))//','//trim(compname(destcomp))//') ', rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if

    end do  ! loop over fields

    if (ESMF_FieldIsCreated(frac_field_dst)) then
       call ESMF_FieldDestroy(frac_field_dst, noGarbage=.true., rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    end if
    if (allocated(data_srctmp)) deallocate(data_srctmp)

    call t_stopf('MED:'//subname)

  end subroutine med_map_FB_Regrid_Norm_All

  !================================================================================

  subroutine map_field_src2dst (fldname, srcfield, dstfield, RouteHandles, mapindex, rc)

    !---------------------------------------------------
    ! map the source field to the destination field
    !---------------------------------------------------

    use ESMF , only : ESMF_LogWrite, ESMF_LOGMSG_INFO, ESMF_SUCCESS
    use ESMF , only : ESMF_LOGMSG_ERROR, ESMF_FAILURE, ESMF_MAXSTR
    use ESMF , only : ESMF_FieldBundle, ESMF_FieldBundleGet
    use ESMF , only : ESMF_Field, ESMF_FieldGet, ESMF_FieldRegrid
    use ESMF , only : ESMF_TERMORDER_SRCSEQ, ESMF_Region_Flag, ESMF_REGION_TOTAL
    use ESMF , only : ESMF_REGION_SELECT
    use ESMF , only : ESMF_RouteHandle, ESMF_RouteHandleIsCreated

    ! input/output variables
    character(len=*)       , intent(in)    :: fldname
    type(ESMF_Field)       , intent(in)    :: srcfield
    type(ESMF_Field)       , intent(inout) :: dstfield
    type(ESMF_RouteHandle) , intent(inout) :: RouteHandles(:)
    integer                , intent(in)    :: mapindex
    integer                , intent(out)   :: rc

    ! local variables
    logical :: checkflag = .false.
    !---------------------------------------------------

    rc = ESMF_SUCCESS

#ifdef DEBUG
    checkflag = .true.
#endif

    if (mapindex == mapnstod_consd) then
       call ESMF_FieldRegrid(srcfield, dstfield, routehandle=RouteHandles(mapnstod), &
            termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       if (dbug_flag > 1) then
          call Field_diagnose(dstfield, fldname, " --> after nstod: ", rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if
       call ESMF_FieldRegrid(srcfield, dstfield, routehandle=RouteHandles(mapconsd), &
            termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_SELECT, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       if (dbug_flag > 1) then
          call Field_diagnose(dstfield, fldname, " --> after consd: ", rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if
    else if (mapindex == mapnstod_consf) then
       call ESMF_FieldRegrid(srcfield, dstfield, routehandle=RouteHandles(mapnstod), &
            termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       if (dbug_flag > 1) then
          call Field_diagnose(dstfield, fldname, " --> after nstod: ", rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if
       call ESMF_FieldRegrid(srcfield, dstfield, routehandle=RouteHandles(mapconsf), &
            termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_SELECT, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       if (dbug_flag > 1) then
          call Field_diagnose(dstfield, fldname, " --> after consf: ", rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if
    else
       call ESMF_FieldRegrid(srcfield, dstfield, routehandle=RouteHandles(mapindex), &
            termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    end if

  end subroutine map_field_src2dst

  !================================================================================

  subroutine norm_field_dest (fldname, dstfield, frac, rc)

    !------------------------------------------------
    ! normalize destination mapped values by the reciprocal of the
    ! mapped fraction or 'one'
    ! ------------------------------------------------

    ! input/output variables
    character(len=*) , intent(in)    :: fldname
    type(ESMF_Field) , intent(inout) :: dstfield
    real(r8)         , intent(in)    :: frac(:)
    integer          , intent(out)   :: rc

    ! local variables
    integer           :: i,n
    integer           :: lrank
    real(R8), pointer :: data1d(:)
    real(R8), pointer :: data2d(:,:)
    integer           :: ungriddedUBound(1)     ! currently the size must equal 1 for rank 2 fields
    integer           :: gridToFieldMap(1)      ! currently the size must equal 1 for rank 2 fields

    ! ------------------------------------------------

    rc = ESMF_SUCCESS

    call ESMF_FieldGet(dstfield, rank=lrank, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    if (lrank == 1) then
       call ESMF_FieldGet(dstfield, farrayPtr=data1d, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       do i= 1,size(data1d)
          if (frac(i) == 0.0_R8) then
             data1d(i) = 0.0_R8
          else
             data1d(i) = data1d(i)/frac(i)
          endif
       enddo
    else if (lrank == 2) then
       call ESMF_FieldGet(dstfield, ungriddedUBound=ungriddedUBound, gridToFieldMap=gridToFieldMap, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       call ESMF_FieldGet(dstfield, farrayPtr=data2d, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
       do n = 1,ungriddedUbound(1)
          if (gridToFieldMap(1) == 1) then
             do i = 1,size(data2d,dim=1)
                if (frac(i) == 0.0_r8) then
                   data2d(i,n) = 0.0_r8
                else
                   data2d(i,n) = data2d(i,n)/frac(i)
                end if
             end do
          else if (gridToFieldMap(1) == 2) then
             do i = 1,size(data2d,dim=2)
                if (frac(i) == 0.0_r8) then
                   data2d(n,i) = 0.0_r8
                else
                   data2d(n,i) = data2d(n,i)/frac(i)
                end if
             end do
          end if
       end do
    end if

    call Field_diagnose(dstfield, fldname, " --> after frac: ", rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

  end subroutine norm_field_dest

  !================================================================================

  subroutine med_map_FB_Regrid_Norm_Frac(fldnames, FBSrc, FBDst, &
       FBFrac, mapnorm, RouteHandle, string, rc)

    ! ----------------------------------------------
    ! Map fldnames in source field bundle with appropriate fraction weighting
    ! ----------------------------------------------

    use ESMF , only: ESMF_FieldBundle, ESMF_FieldBundleIsCreated, ESMF_FieldBundleGet
    use ESMF , only: ESMF_RouteHandle, ESMF_RouteHandleIsCreated

    ! input/output variables
    character(len=*)       , intent(in)           :: fldnames(:)
    type(ESMF_FieldBundle) , intent(inout)        :: FBSrc
    type(ESMF_FieldBundle) , intent(inout)        :: FBDst
    type(ESMF_FieldBundle) , intent(in)           :: FBFrac
    character(len=*)       , intent(in)           :: mapnorm
    type(ESMF_RouteHandle) , intent(inout)        :: RouteHandle
    character(len=*)       , intent(in), optional :: string
    integer                , intent(out)          :: rc

    ! local variables
    integer                        :: i, n
    type(ESMF_FieldBundle)         :: FBSrcTmp        ! temporary
    type(ESMF_FieldBundle)         :: FBNormSrc       ! temporary
    type(ESMF_FieldBundle)         :: FBNormDst       ! temporary
    character(len=CS)              :: lstring
    character(len=CS)              :: csize1, csize2
    real(R8), pointer              :: data_srctmp(:)  ! temporary
    real(R8), pointer              :: data_src(:)     ! temporary
    real(R8), pointer              :: data_dst(:)     ! temporary
    real(R8), pointer              :: data_srcnorm(:) ! temporary
    real(R8), pointer              :: data_dstnorm(:) ! temporary
    real(R8), pointer              :: data_frac(:)    ! temporary
    real(R8), pointer              :: data_norm(:)    ! temporary
    character(len=*), parameter    :: subname='(module_MED_Map:med_map_Regrid_Norm)'
    !-------------------------------------------------------------------------------

    call t_startf('MED:'//subname)
    rc = ESMF_SUCCESS

    call ESMF_LogWrite(subname//' called', ESMF_LOGMSG_INFO)

    call memcheck(subname, 1, mastertask)

    if (present(string)) then
      lstring = trim(string)
    else
      lstring = " "
    endif

    !-------------------------------------------------
    ! Loop over all fields in the source field bundle and map them to
    ! the destination field bundle accordingly
    !-------------------------------------------------

    call FB_reset(FBDst, value=czero, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    do n = 1,size(fldnames)

       ! get pointer to source field data in FBSrc
       call FB_GetFldPtr(FBSrc, trim(fldnames(n)), data_src, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       ! create a new temporary field bundle, FBSrcTmp that will contain field data on the source grid
       if (.not. ESMF_FieldBundleIsCreated(FBSrcTmp)) then
          call FB_init(FBSrcTmp, flds_scalar_name, &
               FBgeom=FBSrc, fieldNameList=(/'data_srctmp'/), name='data_srctmp', rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return

          call FB_GetFldPtr(FBSrcTmp, 'data_srctmp', data_srctmp, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if

       ! create a temporary field bundle that will contain normalization on the source grid
       if (.not. ESMF_FieldBundleIsCreated(FBNormSrc)) then
          call FB_init(FBout=FBNormSrc, flds_scalar_name=flds_scalar_name, &
               FBgeom=FBSrc, fieldNameList=(/trim(mapnorm)/), name='normsrc', rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return

          call FB_GetFldPtr(FBNormSrc, trim(mapnorm), data_srcnorm, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       endif

       call FB_reset(FBNormSrc, value=czero, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       ! create a temporary field bundle that will contain normalization on the destination grid
       if (.not. ESMF_FieldBundleIsCreated(FBNormDst)) then
          call FB_init(FBout=FBNormDst, flds_scalar_name=flds_scalar_name, &
               FBgeom=FBDst, fieldNameList=(/trim(mapnorm)/), name='normdst', rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return

          call FB_GetFldPtr(FBFrac, trim(mapnorm), data_frac, rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       endif

       call FB_reset(FBNormDst, value=czero, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       ! error checks
       if (size(data_srcnorm) /= size(data_frac)) then
          call ESMF_LogWrite(trim(subname)//" fldname= "//trim(fldnames(n))//" mapnorm= "//trim(mapnorm), &
               ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
          write(csize1,'(i8)') size(data_srcnorm)
          write(csize2,'(i8)') size(data_frac)
          call ESMF_LogWrite(trim(subname)//": ERROR data_normsrc size "//trim(csize1)//&
               " and data_frac size "//trim(csize2)//" are inconsistent", &
               ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
          rc = ESMF_FAILURE
          return
       else if (size(data_srcnorm) /= size(data_srctmp)) then
          write(csize1,'(i8)') size(data_srcnorm)
          write(csize2,'(i8)') size(data_srctmp)
          call ESMF_LogWrite(trim(subname)//": ERROR data_srcnorm size "//trim(csize1)//&
               " and data_srctmp size "//trim(csize2)//" are inconsistent", &
               ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
          rc = ESMF_FAILURE
          return
       end if

       ! now fill in the values for data_srcnorm and data_srctmp - these are the two arrays needed for normalization
       ! Note that FBsrcTmp will now have the data_srctmp value
       do i = 1,size(data_frac)
          data_srcnorm(i) = data_frac(i)
          data_srctmp(i)  = data_src(i) * data_frac(i)  ! Multiply initial field by data_frac
       end do

       ! regrid FBSrcTmp to FBDst
       if (trim(fldnames(n)) == trim(flds_scalar_name)) then
          call ESMF_LogWrite(trim(subname)//trim(lstring)//": skip : fld="//trim(fldnames(n)), &
               ESMF_LOGMSG_INFO)
       else
          call FB_FieldRegrid( FBSrcTmp, 'data_srctmp', FBDst, fldnames(n), RouteHandle, rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if

       call FB_FieldRegrid(FBNormSrc, mapnorm, FBNormDst, mapnorm, RouteHandle, rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       ! multiply interpolated field (FBDst) by reciprocal of fraction on destination grid (FBNormDst)
       call FB_GetFldPtr(FBNormDst, trim(mapnorm), data_dstnorm, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       call FB_GetFldPtr(FBDst, trim(fldnames(n)), data_dst, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return

       do i= 1,size(data_dst)
          if (data_dstnorm(i) == 0.0_R8) then
             data_dst(i) = 0.0_R8
          else
             data_dst(i) = data_dst(i)/data_dstnorm(i)
          endif
       end do

       if (dbug_flag > 1) then
          call FB_Field_diagnose(FBDst, fldnames(n), &
               string=trim(subname) //' Mapping (' // trim(fldnames(n)) // trim(lstring), rc=rc)
          if (chkerr(rc,__LINE__,u_FILE_u)) return
       end if

    end do  ! loop over fields

    ! Clean up temporary field bundles
    if (ESMF_FieldBundleIsCreated(FBSrcTmp)) then
       call FB_clean(FBSrcTmp, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    end if
    if (ESMF_FieldBundleIsCreated(FBNormSrc)) then
       call FB_clean(FBNormSrc, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    end if
    if (ESMF_FieldBundleIsCreated(FBNormDst)) then
       call FB_clean(FBNormDst, rc=rc)
       if (chkerr(rc,__LINE__,u_FILE_u)) return
    end if
    call t_stopf('MED:'//subname)

  end subroutine med_map_FB_Regrid_Norm_Frac

  !================================================================================

  subroutine med_map_uv_cart3d(usrc, vsrc, udst, vdst, RouteHandle, rc)

    use ESMF, only : ESMF_Mesh, ESMF_MeshGet, ESMF_MESHLOC_ELEMENT, ESMF_TYPEKIND_R8
    use ESMF, only : ESMF_FieldCreate, ESMF_FieldDestroy, ESMF_FieldRegrid
    use ESMF, only : ESMF_RouteHandle, ESMF_TERMORDER_SRCSEQ, ESMF_REGION_TOTAL

    ! input/output variables
    type(ESMF_Field)       , intent(in)    :: usrc
    type(ESMF_Field)       , intent(in)    :: vsrc
    type(ESMF_Field)       , intent(inout) :: udst
    type(ESMF_Field)       , intent(inout) :: vdst
    type(ESMF_RouteHandle) , intent(inout) :: RouteHandle
    integer                , intent(out)   :: rc

    ! local variables
    integer             :: n
    real(r8)            :: lon,lat
    real(r8)            :: coslon,coslat
    real(r8)            :: sinlon,sinlat
    real(r8)            :: ux,uy,uz
    type(ESMF_Mesh)     :: lmesh_src
    type(ESMF_Mesh)     :: lmesh_dst
    type(ESMF_Field)    :: field3d_src
    type(ESMF_Field)    :: field3d_dst
    real(r8), pointer   :: data_u_src(:)
    real(r8), pointer   :: data_u_dst(:)
    real(r8), pointer   :: data_v_src(:)
    real(r8), pointer   :: data_v_dst(:)
    real(r8), pointer   :: data2d_src(:,:)
    real(r8), pointer   :: data2d_dst(:,:)
    real(r8), pointer   :: ownedElemCoords_src(:)
    real(r8), pointer   :: ownedElemCoords_dst(:)
    integer             :: numOwnedElements
    integer             :: spatialDim
    logical             :: checkflag = .false.
    real(r8), parameter :: deg2rad = shr_const_pi/180.0_R8  ! deg to rads
    character(len=*), parameter :: subname='(module_MED_Map:med_map_uv_cart3d)'
    !-------------------------------------------------------------------------------

    rc = ESMF_SUCCESS

    ! Create two field bundles vec_src3d and vec_dst3d that contain
    ! all three fields with undistributed dimensions for each

    ! Get pointer to input u and v data source field data
    call ESMF_FieldGet(usrc, farrayPtr=data_u_src, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    call ESMF_FieldGet(vsrc, farrayPtr=data_v_src, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! Get pointer to destination data that will be filled in after
    ! rotation back from cart3d
    call ESMF_FieldGet(udst, farrayPtr=data_u_dst, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    call ESMF_FieldGet(vdst, farrayPtr=data_v_dst, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! get source mesh and coordinates
    call ESMF_FieldGet(usrc, mesh=lmesh_src, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    call ESMF_MeshGet(lmesh_src, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    allocate(ownedElemCoords_src(spatialDim*numOwnedElements))
    call ESMF_MeshGet(lmesh_src, ownedElemCoords=ownedElemCoords_src)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! get destination mesh and coordinates
    call ESMF_FieldGet(udst, mesh=lmesh_dst, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    call ESMF_MeshGet(lmesh_dst, spatialDim=spatialDim, numOwnedElements=numOwnedElements, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    allocate(ownedElemCoords_dst(spatialDim*numOwnedElements))
    call ESMF_MeshGet(lmesh_dst, ownedElemCoords=ownedElemCoords_dst)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! create source field and destination fields
    field3d_src = ESMF_FieldCreate(lmesh_src, ESMF_TYPEKIND_R8, name='src3d', &
         ungriddedLbound=(/1/), ungriddedUbound=(/3/), gridToFieldMap=(/2/), &
         meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    field3d_dst = ESMF_FieldCreate(lmesh_dst, ESMF_TYPEKIND_R8, name='dst3d', &
         ungriddedLbound=(/1/), ungriddedUbound=(/3/), gridToFieldMap=(/2/), &
         meshloc=ESMF_MESHLOC_ELEMENT, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! get pointers to source and destination data that will be filled in with rotation to cart3d
    call ESMF_FieldGet(field3d_src, farrayPtr=data2d_src, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    call ESMF_FieldGet(field3d_dst, farrayPtr=data2d_dst, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! Rotate Source data to cart3d
    do n = 1,size(data_u_src)
       lon = ownedElemCoords_src(2*n-1)
       lat = ownedElemCoords_src(2*n)
       sinlon = sin(lon*deg2rad)
       coslon = cos(lon*deg2rad)
       sinlat = sin(lat*deg2rad)
       coslat = cos(lat*deg2rad)
       data2d_src(1,n) = -coslon*sinlat*data_v_src(n) - sinlon*data_u_src(n) ! x
       data2d_src(2,n) = -sinlon*sinlat*data_v_src(n) + coslon*data_u_src(n) ! y
       data2d_src(3,n) =  coslat*data_v_src(n)                               ! z
    enddo

    ! Map all thee vector fields at once from source to destination grid
    call ESMF_FieldRegrid(field3d_src, field3d_dst, RouteHandle, &
         termorderflag=ESMF_TERMORDER_SRCSEQ, checkflag=checkflag, zeroregion=ESMF_REGION_TOTAL, rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

    ! Rotate destination data back from cart3d to original
    do n = 1,size(data_u_dst)
       lon = ownedElemCoords_dst(2*n-1)
       lat = ownedElemCoords_dst(2*n)
       sinlon = sin(lon*deg2rad)
       coslon = cos(lon*deg2rad)
       sinlat = sin(lat*deg2rad)
       coslat = cos(lat*deg2rad)
       ux = data2d_dst(1,n)
       uy = data2d_dst(2,n)
       uz = data2d_dst(3,n)
       data_u_dst(n) = -sinlon*ux + coslon*uy
       data_v_dst(n) = -coslon*sinlat*ux - sinlon*sinlat*uy + coslat*uz
    enddo

    ! Deallocate data
    deallocate(ownedElemCoords_src)
    deallocate(ownedElemCoords_dst)
    call ESMF_FieldDestroy(field3d_src, noGarbage=.true., rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return
    call ESMF_FieldDestroy(field3d_dst, noGarbage=.true., rc=rc)
    if (chkerr(rc,__LINE__,u_FILE_u)) return

  end subroutine med_map_uv_cart3d

end module med_map_mod
